\chapter{EM算法：从原理到代码}

	如果给你一堆数，告诉你这堆数源自某个正态分布，请你利用这堆数来估计这个正态分布的均值和方差。如果这个搞不定，下面的内容估计比较困难，建议先复习极大似然估计，然后再来看接下来的EM算法。
	\section{混合分布的参数估计：一个简单的例子}
\paragraph{动机}	现在有个新的问题，一堆数据，是由两个正态分布混合而成，我们想得到这两个正态分布的参数估计。

\paragraph{建模}为具体地描述这件事，这两个正态分布可以写为，
\begin{align*}
Y_1&\sim N(\mu_1,\sigma_1^2),\;\; \theta_1 = (\mu_1,\sigma_1)\\
Y_2&\sim N(\mu_2,\sigma_2^2),\;\; \theta_2 = (\mu_2,\sigma_2)\\
Y &= (1-\Delta)\cdot Y_1+\Delta\cdot Y_2
\end{align*}

其中，$ \Delta\in \{0,1\} $，且$ Pr(\Delta = 1)=\pi $，即$ \pi $就是数据集中$ Y_2 $存在的概率。这样，整个数据集的密度就可以写成，
\[ g_Y=(1-\pi)\phi_{\theta_1}(y)+\pi\phi_{\theta_2}(y) \]

其中$ \phi_{\theta_1},\phi_{\theta_2} $是不同均值和方差下的正态分布函数密度。那么具有$ N $个样本的对数似然函数就可以写成，
\begin{equation}\label{lh}
 \ell(\theta;data)=\sum_{i=1}^{N}\ln [(1-\pi)\phi_{\theta_1}(y_i)+\pi\phi_{\theta_2}(y_i)]
\end{equation}


其中$ data $表示数据。现在的问题在于上述似然函数很难直接数值最大化，因为对数里面有个和。那么，有这么一个巧妙的办法可以解决它，意即我们要引入一个潜变量。假设对每个观测值$ i $，有一个不可观测的潜变量$ \Delta_i $，它要么取0，要么取1。取0意味着样本来自$ Y_1 $，取1意味着该样本来自$ Y_2 $。那么就可以重新书写似然函数\eqref{lh}式为，
\begin{equation}\label{lh2}
\ell(\theta;data)=\sum_{i=1}^{N} [(1-\Delta_i)\ln\phi_{\theta_1}(y_i)+\Delta_i\ln\phi_{\theta_2}(y_i)]+\sum_{i=1}^{N} [(1-\Delta_i)\ln(1-\pi)+\Delta_i\ln\pi]
\end{equation}

这就把似然函数里面那个和项展开了。但问题在于$ \Delta_i $未知，这个好办，可以用它的期望来代替，即定义，
\[ \gamma_i(\theta)=E(\Delta_i|\theta,data)=Pr(\Delta_i=1|\theta,data) \]

$ \gamma_i(\theta) $通常称为responsibility，$ \theta=(\theta_1,\theta_2)' $。
现在可以按这样的步骤来估计各个参数了。
\paragraph{算法}
\begin{enumerate}
	\item 初始化$ \hat{\mu}_1,\hat{\mu}_2,\hat{\sigma}_1,\hat{\sigma}_2,\hat{\pi} $。
	
	通常情况情况下，随机地取两个观测值$ y_i $作为$ \hat{\mu}_1,\hat{\mu}_2$的初值，然后让$ \hat{\sigma}_1=\hat{\sigma}_2=\sum_{i=1}^N(y_i-\bar{y})^2/N,\hat\pi= 0.5 $。
	\item E步。
	
	利用参数值计算每一个观测值属于某个潜分类，譬如$ Y_2 $时的期望，也即每一个观测值的responsibility，计算公式如下，
	\begin{equation}\label{res}
	 \gamma_i=E(\Delta_i|\theta,data)=Pr(\Delta_i=1|\theta,data)=\frac{\hat{\pi}\phi_{\theta_2}(y_i)}{(1-\hat{\pi})\phi_{\theta_1}(y_i)+\hat{\pi}\phi_{\theta_2}(y_i)}\hspace{2em}i=1,2,\cdots,N
	\end{equation}
	
	为什么这么算？注意到上式是一个基于数据的条件概率。回忆贝叶斯公式，
	\[ Pr(\Delta_i = 1|data,\theta)=\frac{Pr(\Delta_i=1,data|\theta)}{Pr(data|\theta)}= \frac{Pr(\Delta_i=1|\theta)Pr(data|\Delta_i=1,\theta)}{Pr(data|\theta)}\]	
	
	也即通过不需要数据的“无条件概率”$ Pr(\Delta_i=1|\theta) $与特定分类下数据的似然函数$ Pr(data|\Delta_i=1,\theta) $
来搞定。
	\item M步。
	
	用E步得到的$ \gamma_i $，也就是$ \Delta_i $的期望代入那个加权的似然函数\eqref{lh2}式，然后最大化\eqref{lh2}式，可以得到一个$ \hat{\mu}_1,\hat{\mu}_2,\hat{\sigma}_1,\hat{\sigma}_2 $的新的估计，以及把所有点的responsibility进行平均，即$ \hat{\pi}=\sum_{i=1}^N\gamma_i/N $。
	\item E步和M步迭代，直至收敛。
\end{enumerate}

纵观该算法，关键点有这么几个，一是要引入潜变量，算法才好理解。二是要理解responsibility, 即\eqref{res}式。
\section{通常情况下的EM算法}
EM算法主要针对的是直接最大化似然函数很困难，但是通过引入潜变量却可以大大简化计算的情况。

\begin{enumerate}
	\item 初始化参数$ \hat{\theta}_0 $。
	\item 利用初始化的参数计算responsibility。
	\item 根据新的respossibility，重新最大化加权似然函数。
	\item 迭代，直至收敛。
\end{enumerate}

懂了responsibility，至少在程序编写上就没有问题了。
